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We have utilized time-domain magneto-terahertz spectroscopy to investigate the low frequency optical re¬ 
sponse of topological insulator Cu 0 . 02 Bi 2 Se 3 and Bi 2 Se 3 films. With both field and frequency dependence, such 
experiments give sufficient information to measure the mobility and carrier density of multiple conduction chan¬ 
nels simultaneously. We observe sharp cyclotron resonances (CRs) in both materials. The small amount of Cu 
incorporated into the Cuo.o 2 Bi 2 Se 3 induces a true bulk insulator with only a single type of conduction with total 
sheet carrier density ~ 4.9 x 10 12 /cnr and mobility as high as 4000 cnr/V-s. This is consistent with conduction 
from two virtually identical topological surface states (TSSs) on top and bottom of the film with a chemical 
potential ~145 meV above the Dirac point and in the bulk gap. The CR broadens at high fields, an effect that 
we attribute to an electron-phonon interaction. This assignment is supported by an extended Drude model anal¬ 
ysis of the zero field Drude conductance. In contrast, in normal Bi 2 Se 3 films two conduction channels were 
observed and we developed a self-consistent analysis method to distinguish the dominant TSSs and coexisting 
trivial bulk/2DEG states. Our high-resolution Faraday rotation spectroscopy on Cu 0 . 02 Bi 2 Se 3 paves the way 
for the observation of quantized Faraday rotation under experimentally achievable conditions to push chemical 
potential in the lowest Landau Level. 


Topological insulators (TIs) are a newly discovered class of 
materials characterized by an inverted band structure 1 ^caused 
by strong spin-orbit coupling. In the ideal case, they have 
an insulating bulk and only conduct via topologically pro¬ 
tected massless Dirac topological surface states (TSSs). Spin- 
momentum locking in their electronic structure makes TIs 
promising platforms for spintronics application^'. Progress 
in this field has been hampered by the fact that all discov¬ 
ered TIs to date are slightly doped and have a conducting 
bulk. For instance, the proposed topological magneto-electric 
effecP and quantized Faraday rotatiorP remains unobserved. 
Tuning the chemical potential towards the Dirac point and en¬ 
hancing mobility was shown to be very successful in prob¬ 
ing many-body interactions with plasmons and phonons in 
graphene®* and similar advancements are expected in TIs, but 
have not yet been realized. 

The band structure of BLSej is one of the simplest of the 
3D TIs with only a single Dirac cone at the center of the Bril- 
louin zone. Unfortunately, native grown BLSej is known 
to have a conducting bulk due to defects from the growth. 
Suppression of the bulk carrier density has been achieved by 
chemical doping methods 7 ®. Nevertheless, these samples still 
have significant densities of bulk carriers or impurity states 
that are pinned near E/, . Recently it was found that ~ 2% Cu 
doping in thin films suppresses the bulk carriers and allows a 
true insulating state to be realized-. Here we investigate these 
copper doped bulk insulating thin films and their decoupled 
TSSP via magneto-terahertz spectroscopy. These films were 
capped by a thin insulating amorphous Se layer. Details of the 
growth can be found in the supplementary information (SI) 
section P° and in Ref® 


Cyclotron resonance (CR) experiments using THz spec¬ 
troscopy are a powerful tool to study Dirac fermions and probe 
many-body interactions 11 12 . CR is also one of the most accu¬ 
rate measures of effective mass 13 . In previous work, a large 
Kerr rotation in bulk-conducting BLSe 3 films was reported, 
but no obvious resonance was observed^! Cyclotron reso¬ 
nance has been reported in IniSej capped films SI but current 
understanding is that significant Indium diffusion from IniSej 
to BUSejEIIdestroys the simple non-TI/TI boundary at the in¬ 
terface, as a topological phase transition occurs at low Indium 
concentrations (~ 6%R-d 

In the present work, we used time-domain magneto¬ 
terahertz spectroscopy with the polarization modulation 
technique-! (0.5mrad resolution, see SI section I for exper¬ 
imental details 10 ) to observe sharp cyclotron resonances in 
Faraday rotations from both CuoooBLSej and pure BLSej 
thin films. We demonstrate that CuoioBiiSe? can be de¬ 
scribed by a single Drude component with total carrier density 
« 2 d ~ 4.9 x 10 12 /cm 2 . This Drude contribution is consistent 
with pure surface state transport with an E/, ~ 145 meV above 
the Dirac point (75 meV below the conduction band edge), 
which makes Cuo.o 2 BLSe 3 a true topological insulator. CR 
broadening at high field is attributed to an electron phonon in¬ 
teraction. In contrast, we find two channel conduction in pure 
Bi 2 Se 3 films. We determine that the large Faraday rotation 
is induced by a dominant high-mobility TSSs channel with 
an E/. ~ 350 meV. However, a weaker low-mobility second 
Drude term is also required to fit the data. This subdomi¬ 
nant term most likely derives from trivial states (bulk and/or 
2DEG). 

In Fig. |TJa), we compare the zero field THz conductance 
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FIG. 1: (Color online) Data summary for 64 QL Cuo^BFSe, (samplel). (a) Real sheet conductance of 64 QL Cuo.ooBLSej and pure Bi 2 Se 3 
films at 5K. Inset is the Drude conductance Go of 64 QL CuomELSe? after subtracting the phonon and e„ contributions. The arrow indicates 
the deviation from pure Lorentzian form that arises from electron-phonon coupling, (b) Field dependent conductance at 5K. Arrows are guides 
to the eye for cyclotron frequencies, (c) Real and (d) Imaginary parts of the complex Faraday rotation data at different fields at 5K. (e) 
Cyclotron frequency versus field. Solid line is a linear fit. The inset is a cartoon indicating Ep ~ 145 meV above the Dirac point (75 meV 
below conduction band minimum), (f) Complex Faraday angle with fits at 7 T. The solid curve is a fit with a field-dependent scattering rate. 
The dashed curve is a fit using the zero-field scattering rate. Arrows are guides to the eye for cyclotron frequency. 


of a 64 QL Cuo.oaBLSej film to a pure 64 QL BLSe 3 film. 
The Cu incorporated sample’s spectra are characterized by a 
reduced total spectral weight and slightly lower scattering rate 
than the pure BLSe 3 . The spectra can be well fit by an oscil¬ 
lator model with only a Drude term describing free electron¬ 
like motion, a Drude-Lorentz term modeling the phonon and a 
lattice polarizability e m term that originates from absorptions 
above the measured spectral range. 
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Here T’s are scattering rates, u) p ’s are plasma frequencies, and 
d is the film thickness. The spectral weight ( co 2 pD d ) is pro¬ 
portional to the integrated area of each feature in the real part 
of the conductance. It gives the ratio of carrier density to an 
effective transport mass. 
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Here m* is defined as hkp/vp for “massless Dirac fermions”. 
Considering the TSS dispersion up to quadratic corrections 
E = Akp + Bk 2 F , the spectral weight can be expressed in terms 
of kp, where A and B are parameters obtained from ARPES 
(see SI section iP). 


U~pD d = 


k F (A + 2Bkp)e 
2nTi 2 ea 


(3) 


This expression assumes the single channel conduction orig¬ 
inates in two nominally identical TSSs. From the spectral 
weight analysis, having determined k F , we can then calcu¬ 
late both ii 2 D (proportional to k 2 F ) and m*. From these Drude- 
Lorentz fits, we find a total sheet carrier density 112 d ~ 5.0 
xlO I2 /cm 2 , m* ~ 0.14 m e , and E/r ~ 145 (+5) meV according 
to Eqj3] 


We can estimate the carrier density and mass by the zero 
field spectra alone because both can be expressed as a func¬ 
tion of kp. Below we determine the mass in a model free 
fashion through CR experiments. The signature of CR is a 
peak in the real part of the conductance (Fig. |T]b)), an inflec¬ 
tion point in real part of Faraday rotation d' F (Fig. 0c)) and 
a dip in the imaginary part of Faraday rotation 0 F (Fig. 0 d)). 
Its full width at half maximum (FWHM) is the scattering rate. 
Field dependent complex Faraday rotation data are shown in 
Fig. [ljc)(d). One can see that an edge feature around an in¬ 
flection point in 0 F and a dip in 0 F shifts to higher frequency 
with increasing field, which is consistent with CR. We fit the 
data by a Drude-Lorentz model accounting for the field de¬ 
pendence of the Drude term and constraining the parameters 
of the phonon and e a: , by the values extracted from the zero- 
field conductance fits. The conductance in a magnetic field 
can be described by the expression: 
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FIG. 2: (Color online) (a) Real conductance of 100 QL Bi 2 Se 3 comparing single and two component Drude fits as detailed in text. The inset 
enlarges the low frequency regime, (b) Field dependent conductance at 5K. (c) The real and (d) imaginary parts of the complex Faraday 
rotation data at different fields at 5K. (e) The cyclotron resonance frequency versus field. The inset is a cartoon indicating ~ 350 meV 
above the Dirac point, (f) Complex Faraday rotation data with one and two component Drude fits at 7 T. The dashed and solid curves are one 
and two Drude fits respectively. 
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Here the ± sign signifies the response to right/left-hand cir¬ 
cularly polarized light, respectively. oj c is the CR frequency 
to be defined below. The Faraday rotation can be expressed 
as tan(0/r) = -i(t+ - t_)/(t + + f_) 19 . Note that the Faraday 
equation is a complex quantity because, in addition to rota¬ 
tions, phase shifts that are different for right/left-hand polar¬ 
ized light can be accumulated. The imaginary part is related 
to the ellipticitjP 

The fits to this model for the Faraday rotation are shown 
for a representative field of 7 T in Fig. |T|f) (see SI section 
II for fits to all fields 1 ®). In this plot, the dashed curves are 
from a fit with the spectral weight and scattering rate ~ 0.4 
THz set by the zero-field conductance fit (e.g. only a> c al¬ 
lowed to vary). One can see that although the gross features 
of the spectra are reproduced, using the zero field scattering 
rate entirety fails to reproduce certain aspects of the Faraday 
rotation, including the value of the minimum in 0' F . A much 
better fit (solid line) can be obtained by letting the scattering 
rate vary with field, while keeping other parameters (except 
for oj c ) fixed. The origin of this field dependent scattering rate 
will be addressed below. The fits allow us to extract the cy¬ 
clotron frequency as a function of field; it is exhibited by the 
raw spectra as the minimum in 9 F . As shown in Fig. 1 1 le), a 
linear fitting using the expected relation between mass and the 
resonance frequency oj c = eB/m* gives an effective mass of 


0.14 m e . By using Eq. [2] and spectral weight obtained from 
fitting Faraday rotation, we can extract a total sheet carrier 
density « 2 d = 4.9 ± 0.1 x 10 12 /cm 2 . With the known band 
structure of Bi 2 Se 3 , this charge density, and the observation 
of a single kind of charge carrier is only consistent with two 
essentially identical TSSs and an Ep ~145 meV above the 
Dirac point. We can conclude that CuoozBiiSei has an insu¬ 
lating bulk and for the sample highlighted here a high mobil¬ 
ity /j = e/T D m* ~ 4000cm 2 /V- s. The features were robust 
to sample aging as the samples (with the amorphous Se cap) 
maintain bulk-insulating and high mobility after sitting in air 
for eight months (see SI section IP). 

Films of Cuo fliBbSe} can be contrasted with films of pure 
Bi 2 Se 3 , which is known to have the surface chemical potential 
pinned in the bulk conduction band 2021 . We develop a self- 
consistent data analysis method and use our high-resolution 
magneto-terahertz spectroscopy to precisely determine the 
contribution of the subdominant bulk. In Fig. [2ja), we show 
zero field conductance spectra from a typical 100 QL BHSes 
film. It has higher spectral weight than CuoinBFSej, con¬ 
sistent with a higher charge density. For the same reason, 
the plateau-like Faraday rotation at 7 T at low frequencies in 
Fig. [2jf) is as large as -0.25 radians while the value for the 
Cu ( ) Q 2 Bi 2 Se 3 sample is -0.07 radians. 

In Figs. [2jb)-(d), one can see that the CR is exhibited at 
lower frequencies, which indicates that Bi2Se3 has a higher 
CR mass. In Fig. [2je), the linear fit of a> c vs. B gives an effec¬ 
tive mass -0.20 m e . We believe that this derives from TSSs, 
as this mass is inconsistent with the accepted values for the 
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bulk bands 22 or band bending induced surface 2DEG bands 23 . 
Note that the value of CR mass in Bi 2 Se 3 we determine here 
is different than that given in our previous work-31 The reason 
for this discrepancy is discussed at length in the SI section 

VP. 

One can see from Fig. |2jf) that fits of the Faraday rotation 
using only a single Drude term are reasonably good. As before 
we use the spectral weight and CR mass by Eq. [2]to extract a 
total sheet carrier density n^D ~ 1-9 ± 0.1 x 10~/cm 2 . This 
compares favorably to a density of nio ~ 2.0 x 10 13 /cm 2 , a 
mass of ~0.20 m e , and an E/. ~350 meV that we can determine 
purely from an analysis of the spectral weight using the TSS 
dispersion by Eq. [3] We determine a mobility of jj ~ 3200 
cm 2 /V- softheTSSs. 

However, while the fits to the Faraday rotation with a sin¬ 
gle Drude term are excellent, significant discrepancies arise 
if we use the resulting parameters to fit the conductance data. 
In Fig. [2} a), one can see that the single Drude component 
fit significantly underestimates the conductance by a roughly 
constant amount over the entire spectral range. As E/. is in the 
bottom of the conduction band, it is reasonable to ascribe this 
difference to a subdominant low mobility conduction channel 
originating in bulk or 2DEG. Adding a second term with a 
large scattering rate (T/In > 4 THz such that the contribution 
to the real conductance is a constant offset) improves the con¬ 
ductance fit dramatically. If literature values for the bulk or 
2DEG mass are used, this channel has carrier density ri 2 d ~ 
8.0 x 10 12 /cm 2 and low mobility // < 300cm 2 /V-s. The sec¬ 
ond flat Drude term only adds a featureless small background 
to the Faraday rotation. The ratio G\Tsss/Gi to tai > 90% at low 
frequencies confirms that TSSs dominate transport as inferred 
in previous THz measurement®® 7 '. 

Having identified the transport channels in these films and 
their relevant parameters, we can look in more detail at the 
magnetic field dependence of the scattering rates. As shown in 
Fig. [3} a) the scattering rate F/j/2zr in Cuo.o 2 Bi 2 Se 3 increases 
with cyclotron frequency/field, displaying a maximum near 
a cyclotron frequency of 0.9 THz + 0.1 THz. As discussed 
in the SI section III31, we can exclude magnetic field in¬ 
duced spin- or orbital-based electronic mechanisms for this 
broadening®]. We plot the scattering rate vs. the cyclotron 
frequency instead of vs. the magnetic field in order to allow 
a direct comparison with the measurement and phonon fre¬ 
quencies. One can understand magnetic field as a energy sam¬ 
pling process as following. With increasing field, the Drude 
peak center frequency oj c shifts to higher energies. As the 
Drude peak overlaps with the phonon frequency, one expects 
the scattering rate to reach its maximum value and saturate as 
one generally expects that the CR will become broader at fre¬ 
quencies above the relevant phonon frequency because new 
scattering channels emerge. Similar CR broadening has been 
observed in GalnAs quantum wells and graphene and ascribed 
to electron-phonon coupling 25 ®? 

An extended Drude model (EDM) analysis 27 of the zero 
field Drude conductance supports the inference that the broad¬ 
ening is due to coupling to a low energy mode such as a 
phonon. Such analysis allows one to quantify subtle devia¬ 
tions from pure Drude behavior (e.g s. pure Forentzian form) 


( a ) q) c /2ji (THz) 

0.0 0.4 0.8 1.2 



Frequency (THz) 

FIG. 3: (Color online) (a) Scattering rate as a function of cyclotron 
frequency (red). Fully renormalized scattering rate by mass through 
extended Drude analysis as a function of frequency (green). Black 
arrow indicates /f surface phonon frequency at Kohn anomaly, (b) 
Renormalized mass as a function of frequency. The error bars ex¬ 
press the uncertainty in a> pD . 


as shown in Fig. [T] a) inset in the form of frequency dependent 
mass renormalizations m* RN (a>) and scattering rates ro(co)/2n 
shown in Fig(3] In such analysis one inverts the Drude intra¬ 
band contribution to the conductance (Go ) via Eqs. [5] 


r D (co) = 


t(oj) 


^Re(-L), 

4n Go 


mtjco) w; 


m b 
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Here nt/, is the band mass without interaction renormaliza¬ 
tion. An EDM analysis shows a frequency dependent scat¬ 
tering rate To(cu) and an electron-phonon coupling constant 
(T = - l) defined at DC limit ~ 0.55 ± 0.05. Since the 

mass renormalization here is reasonably small, we approxi¬ 
mate the bare plasma frequency with the value extracted from 
zero-field conductance or CR fits. 

From the EDM parameters we can define the “fully renor¬ 
malized” scattering rate T^(w) = which is mani¬ 

fested as the actual width of the low frequency Drude feature 
(e.g. the T D of Eq. [dp? Plotting this quantity in Fig. [3j a), 
we find the general trend and energy scales of the scattering 
rate from the EDM are about the same as those from the CR 
fit. One can see that the result from EDM is a sharper version 
of the curve from the CR fit, which may be expected because 
the CR fit assumed a frequency independent scattering rate, 
and as such means it is an average over a (weakly) frequency 
dependent quantity. 

For the first time we have demonstrated a correspondence 
between the spectral features in a CR experiment and EDM 
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analyses, which indicates the magnetic field itself is not the 
source of the CR broadening. We propose electron-phonon 
coupling as cause of these effects 2 ® ® 132 Jas the energy scale 
of the threshold in the scattering rate matches that of a num¬ 
ber of electron-phonon scales in this system. The value is 
close to the scale of the previously observed Kohn anomaly 
of the surface /j phonon, 0.75 THz at 2k/P®. It also matches 
closely the scale of the maximum acoustic phonon energy that 
can couple to cross Fermi surface scattering c x2kf/27r ~ 0.6 
THz (where c is the acoustic phonon velocity) 31 . The coupling 
strength A extracted by EDM is close to the value A -0.43 by 
analysis of the Kohn anomaly^, and also agrees with calcu¬ 
lation A ~ 0.42 for acoustic phonons in thin film geometrjC 32 . 
The coupling constant we observed is within the different val¬ 
ues ARPES give 33 ®' 7 . The energy scale of the 3 meV mode 
found via ARPES in Ref. 37 ^ agrees well with the characteristic 
energy here. The lower coupling strength we observed may 
be due to the lower chemical potential of our sample than that 
used in Ref.®. 

Looking forward our results demonstrate the possibility of 
observing a quantized Faraday rotation in such films of a 
TI. For the quantized Faraday rotation of a film on a sub¬ 
strate one has tan(0/r) = ®-(v + 1/2) ~ 3.5 mrad(v + 1/2) 


for each surface state in the quantum regime, where a is the 
fine structure constant, v is the filling factor and n is the sub¬ 
strate index of refraction. The Landau level spectrum will 
be E(±v) = +Vf yjlehB | v |. In order to reach the lowest 
LL in this true topological insulator Cuo.oiBESe? will require 
150T, which may be achievable for future THz experiments 
at pulse field facilities. Electric field gating the present sam¬ 
ples in a THz compatible arrangement would bring down this 
field threshold. We hope our high-resolution Faraday rotation 
spectroscopy will enable the observation of the long-sought 
quantized Faraday rotation and topological magneto-electric 
effect®!. 
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Supplementary Information for “High-resolution Faraday rotation and electron phonon coupling in surface states of 
the bulk-insulating topological insulator CuoicBFSej”’ 


I. MATERIALS AND EXPERIMENTAL METHODS 

Standard time-domain THz spectroscopy (TDTS) in a transmission geometry was performed with a custom home-built THz 
spectrometer. In this technique an approximately single-cycle picosecond pulse of light is transmitted through the sample and 
the substrate. The complex transmission is obtained from the ratio of a Fourier transformed time-domain sample pulse over a 
Fourier transformed substrate pulse. The complex conductance can be directly inverted from the transmission equation in the 
thin film limi l 14 ! 17 !: 


f(to) = 


1 + n 

1 + n + ZqG(uj) 


gico(n-\)AL/c 


(i) 


where A L is the small difference in thickness between the sample and reference substrates, n is the real index of refraction of 
substrate and Zq is the vacuum impedance, 376.7 Q. By measuring both the magnitude and phase of the transmission, this 
inversion to conductance is done directly and does not require Kramers-Kronig transformation. TDTS is an ideal tool to study 
the low frequency response of these materials with both the metallic Drude peak and a E|„ infrared active phonon visible. 

Thin films of Cuo.mB^Sej were grown at Rutgers University by molecular beam epitaxy (MBE) on 0.5 mm thick crystalline 
AEO 3 substrates. Details on the growth can be found elsewher^ 2i; . Films grow a quintuple layer (1 unit cell) at a time (1 
QL~ 1 nm). After film growth, a 100 nm amorphous Se cap is deposited to reduce aging effects^. Se capping was shown to 
have a negligible contribution to the optics at THz frequencies®, yet it serves a very important protection layer as Cu doped 
BESe 3 is much less stable in air than pure Bi 2 Se 3 . 3-4% optimal Cu concentration (Cu/Bixl00%) was incorporated during the 
film deposition and the concentration can be controlled at better than the 1% level. Therefore, the ‘x’ in Cu A BiiSe 3 formula 
is 0.015-0.02. Samples were sealed in vacuum and sent immediately to JHU and low temperature TDTS measurements began 
within 24 hours of their growth. The samples were mounted inside an optical helium flow cryostat and cooled down to 5 K 
within an hour. 



FIG. 1: (Color online) Demonstration of experimental setup for Faraday rotation experiments. 

Complex Faraday rotation measurements were performed in a closed-cycle 7 T superconducting magnet at 5 K. We use the 
polarization modulation technique to measure the polarization states accurately®. As shown in Fig. [I] a static wire-grid polarizer 
(WGP1) is placed before the sample. After the sample, a fast rotating polarizer (FRP) unit and another static WGP2 are used. 
WGP1 and WGP2 transmit vertically polarized light. With this combination, in the polarization modulation technique, E x (t) and 
Ey(t) (blue pulses) can be measured simultaneously in a single scan by reading off the in- and out-of-phase outputs from a lockin. 
Complex Faraday rotation Of - 0 F + i0 F can be obtained by (U=atan[£\(«)/F A (w)] after Fourier transforming into the frequency 
domain. Linearly polarized light becomes elliptically polarized after passing through the sample with a Faraday rotation 9 F (real 
part). The imaginary part of the Faraday rotation d ' F is related to the ellipticity in the small rotation angle regime!®. In these 
measurements, a small background rotation from misalignment was subtracted by measuring a blank substrate. 
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II. MORE DATA AND ANALYSIS FOR Cu 0 (cBLSe, AND Bi 2 Se 3 
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FIG. 2: (Color online) Fit quality on 64 QL Cu 0 02 Bi 2 Se 3 (sample 1) at different fields. Arrows are guides to the eye for the cyclotron 
frequencies. Accurate numbers were determined by fits. Solid lines are fits with field-dependent scattering rates. Dashed lines are fits using 
the zero-field scattering rates. 


The carrier density of each TSS can be calculated by the usual relation /i 2 /j = k 2 F /4n. An effective transport mass m* = hkp/vp 
can still be defined for even for ‘massless’ Dirac fermions where the Fermi velocity is determined by v F = dEp/fidk. In our 
analysis we will consider up to quadratic dispersion for surface states (E = Akf + Bk 2 F ) and model the two surface states 
as identical with same carrier density. Considering the TSS dispersion up to quadratic correction, the spectral weight can be 
expressed in terms of k F 39 . 


2 

neo 


f 


Goidu) — — 


kf(A + 2 Bkf)e 2 
2nft 2 eo 


( 2 ) 


Therefore, lower spectral weight means lower kf . Lower k F means lower carrier density and smaller mass. By fitting the zero- 
field conductance, we find (ojj,i,/2Ti) 2 d in the 64 QL Cuo.o 2 Bi 2 Se 3 sample 1 (the one discussed in main text) equal to 3.0+0.2 x 
10 4 THz 2 ■ nm. By using Eq. [ 2 ] we obtain k F = 0.056+0.003 A l , Ep = 145 +5 meV, m* =0.135+0.005 m e and a total sheet carrier 
density ri 2 /j =5.0 ±0.3 xl0 12 /cm 2 . From the Faraday rotation fit, we get the spectral weight (a> p o/2n) 2 d = 2.8 + 0.1 x 10 4 THz 2 - 

nm more accurately. By using the relation co 2 D d = and effective mass from CR experiments, we can get carrier density 
ri 2 D = 4.9 ± 0.1 x 10 12 /cm 2 . These two analysis agree self-consistently. Here we have shown a new way to self-consistently 
study surface state transport. 

Here we show all the data and fits for both 64 QL Cuo.o 2 Bi 2 Se 3 samples (sample 1 and sample 2) . One can see in Fig. [2] 
for Cuo 02 BLSe 3 (sample 1) that the difference in fits between field independent and field dependent scattering rates are highly 
distinguishable at finite fields, which is another indication of increasing scattering rate due to electron-phonon coupling. 

We observed similar phenomena on another 64 QL Cuo.mBLSes (sample 2) and reached a similar conclusion to sample 
1. This sample 2 remained bulk-insulating and with high-mobility after 8 months of exposure to air, which demonstrates the 
realization of a robust bulk-insulating TI by the protection of Se capping. This sample also had ~ 5.0 x 10 l2 /cm 2 but slightly 
lower mobility ~ 2800 cm 2 /V-s. Data and fitting results are shown in Fig. [ 3 ] Here we subtracted the zero-field ‘rotation’ as a 
small background from misalignment, which gives the same result as a substrate background is subtracted. We also performed 
extended Drude analysis and reached a similar conclusion to sample 1. Suppression of the scattering rate at low frequency was 
observed and a coupling constant in the DC limit A ~ 0.50 ± 0.05 is extracted. As discussed in the main text, electron-phonon 
coupling is not sensitive to magnetic field itself. This is what we showed by showing the extended Drude model analysis at 
zero magnetic field gives a overall similar frequency dependence. The magnetic field rather, changes how electron phonon 
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coupling manifests itself in the spectra. One generally expects that such coupling causes broadening of some spectral feature at 
frequencies overlapping with a characteristic phonon energy. By putting on a magnetic field one pushes the zero frequency Drude 
peak to higher frequency and causes it to broaden. The electron-phonon coupling can be seen either by pushing the Drude peak 
across the characteristic phonon frequencies with magnetic field or by looking at subtle deviations from a Lorentzian lineshape 
at frequencies above the phonon frequencies in zero field. Basically, we can view the magnetic field as a energy sampling tool. 




0.0 0.4 0.8 1.2 




Frequency (THz) 



FIG. 3: (Color online) Data summary on 64 QL CuoorBBSes (sample 2) which was exposed in air for 8 months, (a) Complex conductance at 
5 K. (b) Real and (c) Imaginary part of complex Faraday rotation data at different fields at 5 K. (d) Cyclotron frequency versus field. The inset 
is a cartoon indicating E F —150 meV, 70 meV below conduction band minimum, (e) Scattering rate as a function of cyclotron frequency (red). 
Renomalized scattering rate by mass through extended Drude analysis (green) as a function of frequency, (f) Renormalized mass as a function 
of frequency. The error bars express the uncertainty in u> pD . 


In Fig. [4] we show the fits for 100 QL BLSe 3 . This demonstrates that a single channel is the principal contribution to the CR 
and dominates the Faraday rotation. This channel has spectral weight (a> p /2n) 2 d = 7.6 ±0.3 xlO 4 THz 2 ■ nm. We use the spectral 
weight ( ar p d = ^£_) and CR mass 0.20 m e to extract a total sheet carrier density n^o ~ 1.9 ± 0.1 x 10 13 /cm 2 . If using Eq. 
with the TSS dispersion, from spectral weight we find kf ~ 0.11 A -1 , « 2 d ~ 2.0 x 10 13 /cm 2 , m* ~0.20 m e and Ep ~ 350 meV, 
which agrees with ARPES results on the similar samples 4 ^. Agreement between these two analyses self-consistently a ssign this 
channel to topological surface states. Considering that the characteristic penetration depth of the TSS is of order 2 nm 17 ™, ~7 
° rotation per surface is an extremely large rotation per unit length. The second channel gives a small flat background. In the 
CR fits of 100 QL BLSe 3 , we fixed the scattering rate of the low mobility channel to 4 THz as we obtained from zero-field 
conductance fits, but the fits are reasonably insensitive to the precise scattering rate of this channel. Field dependent scattering 
of the TSS channel of the 100 QL BLSe 3 is shown in Fig. [5]and the same basic physics as Cuo.mBLSei was found. Considering 
twice more free parameters in the fits resulted from multiple channel conduction, we are not able to extract field dependent 
scattering of the TSS channel of the 100 QL Bi 2 Se 3 as reliable as the Cuo.mBLSes case. Also because of two Drude coexistence, 
extended Drude analysis is not reliable to extract coupling strength even though suppression of scattering rate at low frequency 
is also observed. Another reason is that it is likely that substantial charge impurity scattering comes from Se defects at the back 
interface and compete with electron-phonon scattering as the former usually give a decreasing scattering rate with magnetic 
field 2 ®. However, electron-phonon coupling is normally expected not to be sensitive to doping and we believe a similar electron- 
phonon coupling exists in undoped BLSe 3 . It is just that the development of truly insulating Cuo.o 2 Bi 2 Se 3 system gives a simpler 
spectra that allows us to determine the field dependent scattering rate systematically and provide quantitative measures of the 
electron-phonon coupling. 

A 32 QL BLSe 3 was also measured. Data summary is shown as Fig. [6] The channel that dominates CR has an effective mass 
m* ~0.19 m e , carrier density U 2 d ~ 1-90 x 10 13 /cm 2 and mobility /r ~2000 cm 2 /V-s, which is consistent with TSSs with Ep ~ 
310 meV above the Dirac point. A second channel with scattering rate ~ 4 THz is added to fit zero-field conductance. By using 
the effective mass of bulk/2DEG from literature, the second channel has total sheet carrier density ~ 0.8 x 10 13 /cm 2 and mobility 
less than 300 cm 2 /V- s. 
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FIG. 4: (Color online) Fit quality on 100 QL Bi 2 Se 3 at different fields. The solid curves are fits using two Drude terms. Dashed curves are 
with one Drude term. Arrows are eye guides for cyclotron frequencies. Accurate numbers were determined by fits. 



FIG. 5: (Color online) Drude scattering rate of 100 QL BLSei as a function of cyclotron frequency at 5 K. 


Note that the Faraday rotation spectra for Bi 2 Se 3 and Cuo.oeBLSei looks similar in shape at first glance, but they are different 
in two ways. Firstly and obviously, CR frequencies is different, which gives a different mass. Secondly and most importantly, 
they indicate two samples are in completely different regime in terms of observation of quantized Faraday/Kerr rotation. Ref. 5 
proposed 90 degree Kerr rotation and Faraday rotation with value of fine structure a ~ 7.3 mrad when chemical potential is 
in the lowest Landau Level (LL) in free standing films. Kerr rotation and Faraday rotation give the same information because 
ultimately it traces back to half-integer quantum hall effect in the conductance. In free standing films in the quantum regime, 
Faraday rotation is tan{6p) = a(v, + 1 /2 + v* + 1 /2), where v, and v/, are LL filling factors for top and bottom surface states. In real 
experiments, free standing films are difficult to reach. People usually measure films on substrate and therefore renormalization 
from substrate effects come in. In the quantum hall regime, taniOic) = ( v / + 1/2 + v* + 1 /2)— 10 mrad x(v, + 1/2 + v* + 1/2), 

taniOf) - + 1/2 + Vb + 1 /2)—3.5 mrad x(v, + 1/2 + v* + 1/2). In this regards. Ref. 14 with 65 degree Kerr rotation has 

never reached the quantum hall regime because chemical potential is way in the conduction band, which is not a topological 
insulator at all. From BQSea to Cuo.ccBLSes, we pushed LL filling factor v ~ 40 to v <10 with 7 Tesla magnet. The Laudau 
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FIG. 6: (Color online) Data summary on 32 QL Bi 2 Se 3 at 5 K. (a) Complex conductance at 5 K. Solid lines are 1 Drude fits using parameters 
obtained from fitting the Faradary rotation with 1 Drude term. Dashed lines are 2 Drude fit. (b)Real and (c) Imaginary Faraday rotation at 
different fields, (d) Cyclotron frequencies as a function of field. Solid line is a linear fit of u> c = eB/m*. The inset is a cartoon indicating 
E f —310 meV. 


Level energy spectrum is E(±v) = ±\’p eTiB \ v |, where ± stands for electron and holes representatively. With experimentally 
achievable magnet field such as pulse magnet facility at Los Alamos National Laboratory, 150 T is enough to push the chemical 
potential into lowest LL and observe the quantized Faraday rotation. We just need to wait the THz system with pulse field 
finish its construction. Recall that scattering rate in CuQQ^BLSet is only less than 2 meV, so LLs will be well separated in 
energy. The threshold for magnetic field to reach lowest LL can be reduced further by combining with gating. Also, we 
reduced resolution from 5 degree in Ref.® to 0.03 degree (0.5 mrad) in the current work. As discussed in the main text, our 
high-resolution Faraday rotation spectroscopy on Cuo.mBLSea opens the new era to observed quantized Faraday rotation and 
topological magneto-electric effect, which are smoking guns for Dirac surface states. 


III. TRANSPORT LIFETIME OF TOPOLOGICAL INSULATORS IN AN EXTERNAL MAGNETIC FIELD 


In this section we investigate various possibilities for the increase in scattering rate at high magnetic fields deriving from 
effects due to the field itself. 


A. Zeeman (Spin) Effect On Transport Lifetime Of Topological Insulators In An External Magnetic Field 

In principle the canting of spins in the surfaces states due to Zeeman coupling under applied magnetic field can cause an 
increase in the scattering rate because back scattering can occur. In practice this is a very small effect at the current chemical 
potential levels. One can see this as follows. 

In the topological regime where the bulk is insulating, one can describe the light-matter interaction of a topological insulator 
by the Hamiltonian of its surface states: 


H k = vcr • (k x z), (3) 

where k is the electron momentum, v is the Fermi velocity, and cr is the Pauli matrix corresponding to the spin degrees of 
freedom. In an external out-of-plane magnetic field B = V x A (where A the corresponding vector potential), the Hamiltonian 
Eq. (|3]> couples to the magnetic field through the orbital degrees of freedom with k —> k-eA/c , and the spin degrees of freedom 
with H k —> //ft + A where the Zeeman coupling A = -g s /j B B/2 acts like a Dirac mass term. We examine the effect due to each 
coupling on the transport lifetime. 
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In low magnetic fields, the effect of magnetic field on the transport lifetime is largely due to the spin polarization effect from 
Zeeman coupling. Therefore we first consider the regime a> c T tr , [A|r„ <s 1, where u> c = eBv/Hk F is the cyclotron frequency. It 
will be convenient to derive the Zeeman field dependence of the transport lifetime T tr using the quantum kinetic theorjEl. The 
kinetic equation for the momentum-dependent charge-spin density matrix f k is 


dfk 

dt 


+ i[H k + Act ,, f k ] + 


1 

2 


e 

eE + -v k x B 

c 





(4) 


where E is the external A.C. electric field from the incident light, v k = dH k /dk = v(zxcr) is the single-electron velocity operator, 
and [ , ] + is the anti-commutator. J(f k \k,t) denotes the collision integral containing information about scattering rates. Note 
that Zeeman coupling mixes different spins and gives rise to a “spin canting” effect, opening up a spin-flipping channel through 
electron-impurity scattering. A related problem was studied in the literature in the context of the anomalous Hall effect of 
topological insulators due to magnetic impurities®. 

We first start by decomposing the density matrix as f k = ff }> + 6f k , where / A <0) is the equilibrium distribution function and 
6fk is the non-equilibrium part. We can facilitate the solution of Eq. (|4ji by resolving 6f k into components of projected onto the 
charge and spin sectors: 6f k = n k + s k ■ cr, where n k and s k describe the charge and spin density distributions respectively. For 
the semiclassical weak-field regime we are studying in this section, the Drude response is given by the conductivity to leading 
order in 1 /E F T tr , and the transport lifetime is determined only by the dynamics of n k and the component of s k projected along 
the electron’s spin direction at the momentum k. We find the momentum relaxation rate as 


1 

2T k 


ths k r 

4v 2 Jo 


—— \u(k, <p k ’ k )\ (1 - cos (f> k ' k ) ^1 + cos- 0 k + sin" Q k cos 



(5) 


where ipk'k is the angle between the momenta k and k ', sin 9 k = vk/e k and cos 6 k = A/e k with s k - ^(vk ) 2 + A 2 , «,■ is the impurity 
concentration, and \u(k, (p k - k )\ 2 is the disorder-averaged impurity potential evaluated at momentum k' = k. In the absence of a 
magnetic field such that A = 0, the angular form factors in Eq. {5]) reduce to the standard expression (1 - cos 2 <p k ' k ) which 
describes the suppression of back scattering of TSS. 

We can further model the impurity potential as short-range delta-function scatters with | u(k, (f> k ’ k )\ 2 = u 2 being independent 
of k. The transport lifetime is determined by Eq. ([5]» at the Fermi level at low temperatures k F T <sc E F . Evaluating the angular 
integral, we then obtain the transport lifetime in the weak magnetic field regime 


1 

2r fr 


njU^Ef 


8 v 2 


1+3 


8 s BbB 
2 E F 


( 6 ) 


Denoting the zero-field transport scattering rate as 1 /2 t„ q = n,n 2 E F /ttv 2 , the weak-field transport scattering rate 1/2 T tr ,B is 
related to that in zero field as 


1 

Ttr.B 


1 

Tir.O 


1+3 


8sBbB \ 2 

2 E F ) 


(7) 


If we use g, ~ 50, E F ~ 150meV, the spin effect only gives a 0.27% increase at 3 T, which is orders of magnitudes off from 
experiments. Therefore, we do not believe the Zeeman effect is the cause for the increase of scattering rate. 


B. Orbital Effect On Transport Lifetime Of Topological Insulators In An External Magnetic Field 

In this section, we consider the effect of orbital coupling to the magnetic field on the transport lifetime. The influence of 
impurity scattering on the density of states and magneto-transport properties of the electron system can be captured using the 
self-consistent Born approximation (SCBAj®. This approximation captures certain features of Landau level (LL) broadening 
due to disorder and is often sufficient to describe magnetotransport properties when localization effects are not important, as 
expected in optical experiments. 

The SCBA for the case of graphene was discussed extensively in Ref.®. For our Cu-doped sample, we estimate from Faraday 
rotation and SdH oscillation^ that the minimum occupied LL which occurs at the largest field li - 7 T is n ~ 10. Therefore, 
the Fermi level is far from the 11 — 0 Landau level (LL), and the effect of inter-LL coupling between opposite levels n and -n 
is negligible. The LL spacing at the Fermi level is of the order of meV, which coincides with the energy scale of the scattering 
rate we extracted from our data [see Fig. 3(a) in main text]. This seems to suggest the sample would be in the regime where 
o.) r T q < 1, where l/r q is the single-particle lifetime that captures the LL disorder broadening (to be distinguished from 1 /r„. 
which is the transport lifetime we want to calculate. Generally l/r ? > 1 /r, r .) 
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The self-energy matrix is diagonal in the LL index and given by 


N c 

2 Y 7 




e-E 




( 8 ) 


where we have defined the dimensionless parameter p = n,^/(Am-' 1 ) that characterizes the disorder strength. Eq. (|8| defines the 
self-consistency condition from which the self-energy can be solved. For high LL filling |A| » 1, one can use the approximation 
of extending the lower limit of the sum to —N c and setting N c —> oo, because the main contribution comes from the Fermi level 
and contributions from negative LLs and from N to oo amount to negligible errors. 

Performing the sum in Eq. ([8]» using the Poisson summation formula, 


E = np (e - E) cot | -j [(e - E) 2 - A 2 ] }, 
and then expanding the cot function into a Fourier series, we obtain 


k= 1 


E' + iE" = -inp (e - E' - iE") <{ 1 + 2 ^ 4 ex P \'[( e “ s ') 2 - A 2 - (E") 2 ] [ [, 


where we have denoted the real and imaginary parts of the self-energy as E = E' + /E". A is the so-called Dingle factor 

A d = exp tIE"! 0 - E')l. 


(9) 


( 10 ) 


( 11 ) 


Eq. p| is a nonlinear equation that needs to be solved iteratively. Ad serves as the small parameter with which Eq. ( | 1 ()[ > can be 
solved for E up to first order in disorder correlation ~ n,u 1 )) and to first order in Ad (0(A 1 d )). 

First, in the zeroth order 0(A° D ) we have E' = 0 and E" = -(?ijU 2 J4v 2 )e. The latter can be identified as the quasiparticle 
scattering rate 1/2 t„q = (n,nJ/4v 2 )e at zero magnetic field. 


To first order 0(A x n ), we obtain 


E' = -eTosin 

T ?,o 


2 n 


e 2 - A 2 


E" = -- 


2 r, 


q, 0 


1 + 2A d cos 


2n 


e 2 - A- 


The density of states is given by v(s) = 21m l./a> z p, therefore we obtain from Eq. (12 1 


v(s) = vq( s) < 1 + 2A d cos 


2 n 


s 1 - A 2 


( 12 ) 


(13) 


with vo(s) = s/2nv 2 being the density of states at zero magnetic field. 

Since the ratio of the transport scattering rate at finite field to that at zero field is equal to that of the density of states 4 -^, finally 
we obtain 


1 _ 1 

Ttr.B T tr Q 


1 + 2 A d COS 


2n 


E\ - (g s pBB/2) 2 


(14) 


where l/r fr o = H,My£’/r/4v 2 is the transport time at zero field. 

In this scenario, the scattering rate should show magneto-oscillations above 3 T when a> c Ttr,B > 1. As we saw no such 
oscillations and only an increase in the scattering rate with field, we also rule out this possibility. 


IV. EFFECT OF Se CAPPING 

Our previous work on pure BESes showed that Se capping had a negligible effect on the optical properties of BESe), but 
does decrease the scattering rate somewhat 17 . However, it is of great importance in the present case of Cu doped Bi 2 Se 3 as these 
samples are very air sensitive. Transport measurements in Ref. 9 were performed right after samples were taken out of ultra-high 
vacuum MBE chamber. Samples for THz measurements were sealed in vacuum bags and shipped overnight to JHU. The total 
exposure to atmosphere is around 30 mins before loading into the cryostat and plus the overnight in the raw-vacuum bag. As 
one can see in Fig. [TJa), Se capping protects the mobility of the sample well as opposed to non-capped samples with a flat Drude 
component. lOOnm Se capping still protects the sample well even after 8 months explosure in air as shown in Fig|7|b) on sample 
2 . 
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FIG. 7: (Color online) (a) Real conductance of 64 QL Cuo 02 BbSe 3 (sample 2) with and without lOOnm Se capping at 5 K . (b) Real conductance 
of 64 QL Cu 0 . 02 Bi 2 Se 3 (sample 2) with lOOnm Se capping measured right after receiving the sample and 8 months later. 


V. LIMITATIONS OF USING ONE DRUDE COMPONENT FIT IN PREVIOUS WORkPHI 

For the 100 QL BLSe 3 sample, if we fit the zero field conductance data by only a single Drude term, a phonon term and a e m 
term, we find the Drude spectral weight (a> p p,l2T:) 2 d is 1.25 x 10 5 THz 2 ■ nm. Using Eq. § one gets kp ~0.14 A *, m* ~0.22 
m e and E h ~480 meV. We believe these values are overestimated and the low frequency spectral weight is actually comprised 
of two terms. As detailed above, a fit to the Faraday rotation data by one Drude term alone shows that the spectral weight 
((LtpD/2n) 2 d contributing to the TSSs is only 7.6 x 10 4 THz 2 -nm. Using Eq. [ 2 ] kp —0.11 A -1 , m* -0.20 m e and Ep ~ 350 meV 
are obtained, as already mentioned above. One can see Ep and k F are overestimated by ~ 30%, while in* is overestimated by ~ 
10% if one associates all the Drude spectral weight in the zero-field conductance to the TSSs. These quantities are overestimated 
in the single Drude component model as one assigns all the bulk/2DEG spectral weight as the TSSs spectral weight. The second 
Drude channel (bulk/2DEG) is flat in our measurement frequency regime and only contributes a smooth background. However, 
considering the fact that the TSSs contribute more than 90% to the total conductance at low frequencies, a single Drude model 
analysis is still a good approximation and a thickness independent Drude peak as observed previously in BLSeJ® 17 is not 
surprising. The presence of this very flat subdominant contribution does not call into question any of the conclusions of our 
previous works^ 17 . Here we just clarify this point and provide a method to separate contributions from TSSs and bulk/2DEG 
in magneto-THz measurements. 


VI. EFFECTIVE MASS INCONSISTENCY WITH PREVIOUS WORKP 

In our earlier work 1 “H based on a fit with just a single Drude term, a heavier 0.35 m e cyclotron mass for TSSs was reported. 
This was judged to be a reasonable number based on a linear dispersion of the TSSs as it gave a reasonable Ep ~0.5 eV using 
m* = Ep/v 2 F and v F = 5 x 10 5 m/s. However, there are a number of inconsistencies that this mass presents. If one considers 
quadratic corrections to the TSS dispersion, one finds that this mass does not give a consistent kp. Moreover, if one uses a linear 
dispersion Ep = Tivpkp and uses the spectral weight of a single Drude component to estimate the chemical potential, the relation 
is: 


2 , Tivpkpe 2 E F e 2 


27T€o Ti 2 2neoti 2 


(15) 


Here we assume two almost identical surface states. An analysis of the spectral weight gives Ep ~0.75 eV, which is unreasonably 
high. Therefore, when Ep is in the conduction band, quadratic corrections must be significant. One may note that in order that 
the wave functions of the TSSs vanish at the Tl/vacuum interface, a quadratic correction is required 4 ®!. | n contrast, when we 


use Eq. 15 based on linear dispersion to estimate the chemical potential of Cuo.mBQSes, Ep is found to be 170 meV + 10 meV. 
This better agreement is expected because surface state dispersions deviate from linearity less when Ep is closer to Dirac point. 
In summary, we believe that the CR mass 0.35 m e given in previous work Ref.® was incorrect for pure BioSes films and the 
revised value of 0.19-0.20 m e given in this work is correct. 
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FIG. 8: (Color online) Calculation of various measures of the complex Kerr rotation with effective mass 0.35 m e and mass 0.20 m e with carrier 
density n 2 o = 2 x 10 l3 /cm 2 and scattering rate f D =l THz at 6 T. 


The reasons for this correction is multifold. 1.) Signal to noise has been improved by orders of magnitude in the current 
generation of experiments with polarization modulation technique with 0.5 mead resolution-^ allowing more reliable fits to the 
spectra. 2.) The current generation of BtySe 3 films have mobilities of the order 2000 - 3200 cm 2 /V- P. In order to see a 
well-defined cyclotron resonance in the dissipative response, the relation lj c t = iiB > 1 needs to be satisfied. oj, t = juB ~ 1.5-2 
is still too small to see very sharp peaks in the amplitude of Kerr/Faraday rotation in our field range. 3.) In the previous work, 
due to (now overcome) uncertainties in the measured phase in the polarization modulation experiment, only the amplitude of 
the Faraday/Kerr rotation was considered reliable. Unfortunately this is a somewhat insensitive method of probing the cyclotron 
frequency. As one can see from Fig. [8] the rotation spectrum does not differ much using 0.35 m e or 0.20;;;,, in the amplitude 
plot. To fit the data accurately when oj c t ~ 1, one needs to fit the complex Faraday rotation; the peak in imaginary part (-0" K ) 
approximately gives the direction evidence for cyclotron resonance. RefPdid not observe obvious signature for CRs because 
of measuring amplitude of Kerr rotation with poor resolution. In this regard, phase sensitive time-domain THz spectroscopy is 
a powerful tool to study cyclotron resonances in topological insulators of current generation. Despite the differences in these 
numbers, the Bi 2 Se 3 samples we measured in this paper are similar to those in RefBJ and the qualitative conclusion that TSSs 
dominates made in ReP remains. 


VII. ASSUMPTION OF NEARLY EQUAL CONTRIBUTION OF T HE TOP AND BOTTOM SURFACE STATES 


In the main text, we showed that the carrier density estimated by using the spectral weight and mass from the Faraday rotation 
fits is close to the value obtained using spectral weight from Faraday rotation fit and surface state dispersion. This shows the 
assumption of nearly equal contributions from the top and bottom surface states is a good approximation. Similarly, DC transport 
have both shown the chemical potential of the two surface states differed by around 50 meV^S. 

A method was discussed in Ref. 15 to distinguish the contribution of the two surfaces from each other. Unfortunately, it appears 
that the interface between In 2 Se 3 and Bi 2 Se 3 which present in the In 2 Se 3 capped samples used in that work is not as simple as 
assumed in Ref. 15 . It was been found that there is 20-30% Indium diffusion into the BBSes layer in a recent stud\P The 
topological phase transition occurs near x~ 6 % in (Bi ^Irn^SeJ 17 51 • Therefore, the interface of In 2 Se 3 and BBSes is not 
the boundary of normal band insulator and topological insulator. Indium probably diffuses into BBSes layer with a gradient. 
The true interface TSS must be buried deeply and exist in a background of high disorder. Due to Indium diffusion, the interface 
which hosts SSs could be (Bio.o 6 lno. 94 ) 2 Se 3 / (Bio.o 4 sIno. 955 ) 2 Se 3 . This hypothesis may explain that the conduction band minimum 
is positioned only ~ 80 meV above the Dirac point, because the bulk gap decreases when approaching the topological phase 
transition point at x ~ 6 %. Also, more Indium substitution reduces the total carrier density 51 , which could be the reason that the 
top surface has a lower carrier density as discussed in Ref. 15 . We also measured a 32 QL Bi 2 Se 3 film capped by 10 nm In 2 Se 3 . 
This sample has lower spectral weight and a larger scattering rate, as shown in Fig. [9] 
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FIG. 9: (Color online) Real conductance of 32 QL BF>Se 3 with and without lOnrn ImSe^ capping at 5 K. 


VIII. ROLE OF COPPER 

The majority of Cu was found to be Cu 0+ (neutral) by X-ray photo-emission spectroscopy in this batch of sampled Our 
data supports this scenario. If 2% Cu substitutes for Bi, then we would observe the shift of phonon frequency as we did in the 
(Bii„ v In v ) 2 Se 3 case. In (Bi] v In v ) 2 Se 3 , we can observe a shift of the phonon at x~ 1%. Also note that Cu has around half of the 
atomic number of In, which means 2% Cu substitution would shift the phonon frequency more than In substituted case. Normal 
Bi 2 Se 3 is known to have a conducting bulk with chemical potential 350-450 meV due to Se vacancies. In keeping with RefP, 
we believe Cu incorporation passivates Se vacancies, so the carrier density is reduced. We encourage more work to be done to 
resolve the role of copper by other methods such as STM measurements and first-principle DFT calculations. 


IX. BAND BENDING 

Topological Insulators can be understood as special narrow-gap semiconductors with surface states. Therefore, band bending 
effects can be important in such systems. Ref . 9 has a nice summary about band bending in TIs. Here we re-emphasize its 
importance. The bulk chemical potential for normal BESe 3 is pinned near the conduction band minimum (therefore ~ 220 
meV above Dirac point while the surface state chemical potential is ~ 350 meV in our samples. Downward band bending 
results in accumulation layers. From magneto-THz measurements, we concluded that any accumulation layer carriers have a 
large scattering rate, are essentially featureless in the Faraday rotation and count for less than 10 % of total conductance at low 
frequencies. For Cu 0 . 02 Bi 2 .Se 3 , the chemical potential at the surface is ~ 150 meV above the Dirac point, and therefore upwards 
band bending must occur. ARPES can play an essential role in determining whether a material is a topological insulator by 
counting if the number of surface state branches is odd or even, but it is not the best tool to conclude whether a TI is bulk- 
insulating or not due to its extreme surface sensitivity and band bending. For example, upwards band bending was reported 
in bulk-conducting TIsP where ARPES did not observe bulk states but SdH saw a dominating bulk contribution. The most 
effective way to probe bulk-insulating TIs is through transport measurements. In DC transport, if the carrier density contributing 
to SdH oscillations is equal to the total carrier density measured by the Hall effect, then the bulk is insulating. In AC optics, 
if the carrier density contributing to cyclotron resonance is equal to the total carrier density in the Drude term of the zero field 
conductance, one can conclude that the bulk is insulating. 





